* VAR code for results in Leduc, S. and K. Sill, "Expectations and Economic Fluctuations: An Analysis Using Survey Data,"
* Review of Economic and Statistics 95(4), October 2013, 1352-1367.
* This code generates the responses to a positive shock to expected unemployment in in the first and third columns of Figure 8, which uses the Livingston expectation data.
* To reproduce the results in Figure 8, the positive expectation shock must be converted to a negative one.
* The different expectations in columns 1 and 3 can be changed on line 257. The estimation period also needs to be changed on lines 272-273.


********************************************************************************
*                                                                              *
* Note:  Computes 3 Confidence Intervals Simultaneously (90%, 80%, 68%)        *
*                                                                              *
********************************************************************************

********************************************************************************
* Note:  This program renormalizes the impulse responses to VAR shocks
*        by dividing all impulse responses by the own-contemporaneous response
*        (see the two "Patch Sections")
*
*        This program also computes historical decompositions & structural shock
*        estimates.
********************************************************************************

*********************************************************************************************
* Program Purpose:
* (1) Reads in semiannual Livingston survey data:
*     yy:1 = the June survey (which forecasts DEC)
*     yy:2 = the Dec  survey (which forecasts JUNE)
*
* (2) Reads in monthly macro variables as semiannual data and constructs all transformations:
*     yy:1 = avg of May(5,t), June(6,t), July(1,t+1), Aug(2,t+1), Sep(3,t+1), and Oct(4,t+1)
*     yy:2 = avg of Nov(5,t), Dec(6,t),  Jan(1,t+1),  Feb(2,t+1), Mar(3,t+1), and Apr(4,t+1)
*
* (3) VAR: pdot, expec, u, r
*     IDENTIFICATION:  Cholesky schemes
*
* (4) Computes optimal lag selection via AIC, SIC, but also allows for arbitrary lag choice
*
* (5) Estimates VAR
*
* (6) Computes impulse responses, bootstrap confidence intervals for VAR shocks, with all
*     impulse responses renormalized by dividing by the contemporaneous own response
*
*********************************************************************************************

cal 1947 1 2
alloc 0 2011:1


****************************************
*                                      *
* Part 1. Read Livingston Survey Data  *
*                                      *
****************************************

open data c:\myfiles\papers\news\code\data_restat.rat
data(format=rats,org=obs) 1947:1 2011:3 ruc12med rgdp10yrmed

****************************************************************
*                                                              *
* Part 2. Input Macro Variables                                *
*         As Semi-Annual Observations                          *
*                                                              *
****************************************************************


*******************************
* 1. Jan & July Monthly Values
*******************************
data(format=rats,org=obs,select=1) 1947:1 2011:3  ruc  rtb3 cpi
set ruc1      / = ruc(t)
set r1        / = rtb3(t)
set p1        / = cpi(t)

clear ruc rtb3 cpi

*******************************
* 2. Feb & Aug Monthly Values
*******************************
data(format=rats,org=obs,select=2) 1947:1 2011:3  ruc  rtb3 cpi
set ruc2      / = ruc(t)
set r2        / = rtb3(t)
set p2        / = cpi(t)

clear ruc rtb3 cpi

*******************************
* 3. Mar & Sept Monthly Values
*******************************
data(format=rats,org=obs,select=3) 1947:1 2011:3  ruc rtb3 cpi
set ruc3      / = ruc(t)
set r3        / = rtb3(t)
set p3        / = cpi(t)

clear ruc rtb3 cpi

*******************************
* 4. Apr & Oct Monthly Values
*******************************
data(format=rats,org=obs,select=4) 1947:1 2011:3  ruc  rtb3 cpi
set ruc4      / = ruc(t)
set r4        / = rtb3(t)
set p4        / = cpi(t)

clear ruc rtb3 cpi

*******************************
* 5. May & Nov Monthly Values
*******************************
data(format=rats,org=obs,select=5) 1947:1 2011:3  ruc  rtb3 cpi
set ruc5      / = ruc(t)
set r5        / = rtb3(t)
set p5        / = cpi(t)

clear ruc rtb3 cpi

*******************************
* 6. Jun & Dec Monthly Values
*******************************
data(format=rats,org=obs,select=6) 1947:1 2011:3  ruc  rtb3 cpi
set ruc6      / = ruc(t)
set r6        / = rtb3(t)
set p6        / = cpi(t)

clear ruc rtb3 cpi

************************************
* 7. Reading data for real-time RGDP
************************************

*******************************
* First and third quarters
*******************************
data(format=rats,org=obs,select=1) 1947:1 2011:3 rgdprt2000q2 rgdprt2007q2 rgdpunrevised grgdpunrevised grgdp
set rt2000rgdp1  / = rgdprt2000q2(t)
set rt2007rgdp1   / = rgdprt2007q2(t)
set rgdpu1       / = rgdpunrevised(t)
set grgdpu1       / = grgdpunrevised(t)
set grgdpr1         / = grgdp(t)

clear rgdprt2000q2 rgdprt2007q2 rgdpunrevised grgdpunrevised grgdp

*******************************
* Second and fourth quarters
*******************************

*data(format=rats,org=obs,select=2) 1948:1 2007:2 rgdprt2000q2 rgdprt2007q2 rgdpunrevised grgdpunrevised
data(format=rats,org=obs,select=2) 1947:1 2011:3 rgdprt2000q2 rgdprt2007q2 rgdpunrevised grgdpunrevised grgdp
set rt2000rgdp2  / = rgdprt2000q2(t)
set rt2007rgdp2   / = rgdprt2007q2(t)
set rgdpu2       / = rgdpunrevised(t)
set grgdpu2       / = grgdpunrevised(t)
set grgdpr2       / = grgdp(t)

clear rgdprt2000q2 rgdprt2007q2 rgdpunrevised grgdpunrevised grgdp

*************************************************************************************
*                                                                                   *
* Part 3. Compute VAR Semi-Annual Macro Data                                        *
*                                                                                   *
* Note 1. The first period  in each year (yyyy:1) is the Livingston June Survey     *
*         The second period "   "     "  (yyyy:2)  "  "     "       Dec  Survey     *
*                                                                                   *
*************************************************************************************

* Unemployment Rates
set urate   /  = (ruc5(t)  + ruc6(t)  + ruc1(t+1)  + ruc2(t+1)  + ruc3(t+1)  + ruc4(t+1))/6.0

* Interest Rates
set r       /  = (r5(t)    + r6(t)    + r1(t+1)    + r2(t+1)    + r3(t+1)    + r4(t+1))/6.0

* Inflation Rates
set inf     /  = (log(p4(t+1)/p4(t)))*2.0*100.0


* Real-time real gdp
set rgdprt    / = log((rgdpu2(t) + rgdpu1(t+1))/2.0)

* Real-time real gdp growth
set grgdprt    / = ((grgdpu2(t) + grgdpu1(t+1))/2.0)

* Revised real gdp growth
set grgdprev    / = ((grgdpr2(t) + grgdpr1(t+1))/2.0)


*Oil shock quantitative dummy
set oshock  /  = 0.0
set oshock 1973:2 1973:2 =  7.8
set oshock 1978:2 1978:2 =  8.9
set oshock 1980:1 1980:1 =  7.2
set oshock 1990:1 1990:1 =  8.8

set oshockreplace  / = oshock

*Fiscal policy dummy
set gshock  /  = 0.0
set gshock 1965:1 1965:1 = 1
set gshock 1980:1 1980:1 = 1
set gshock 2001:2 2001:2 = 1

clear  ruc1     ruc2     ruc3    ruc4    ruc5    ruc6
clear  r1       r2       r3      r4      r5      r6
clear  p4       pcore4
clear rt2000ip1 rt2000ip2 rt2000ip3 rt2000ip4 rt2000ip5 rt2000ip6
clear  ism1     ism2     ism3    ism4    ism5    ism6
clear rt2000rgdp1 rt2000rgdp2 rt2007rgdp1 rt2007rgdp2

************************************************************
*                                                          *
* Part 4.  Program Parameters For Lag Selection Routine    *
*                                                          *
************************************************************



*******************************************
* Set Generic Variables to Use in the VAR
*******************************************

set pdot      /  = inf;         * actual inflation measure
*set expec     /  = ruc12med;   * expected unemployment rate
set expec     /  = rgdp10yrmed; * expected unemployment rate
set u         /  = urate;       * unemployment rate measure
set r         /  = r;           * nominal rate measure
set y         /  = grgdprt;     * real-time output indicator

*********************************************************************************************
* Set Sample Periods & Other Parameters
*
* Note 1. We do not want lags of RHS variables to extend into the period before 1952:1
*         (in the first sample) or into the the period before 1979:2 (in the second sample).
*         Thus, given the samples listed below, we adjust the starting point of each sample
*         in the estimation commands, depending on the number of lags used in the VAR.
*********************************************************************************************

*compute sbeg1 =  1965:2;
*compute send1 =  2008:2;

compute sbeg1 =  1991:1;
compute send1 =  2010:1;

declare rectangular[integer]  lagchoice(1,3)

        * Column 1  :   Corresponds to AIC (determined below)
        * Column 2  :   Corresponds to SIC (determined below)
        * Column 3  :   Any number you want

compute lagchoice(1,3) =  2 ; * enter an integer for number of lags to use, for use when you do
                              * not want to use the AIC or SIC (pre-Volcker period)

compute LS1   = 3;          ; * Lag Scheme(1 = AIC, 2 = SIC, 3 = other)

compute nvar  = 5           ; * number of variables in the VAR

                              * List order of variables as given below in the VARIABLES command

compute [vector[string]] varnames = ||'pdot', 'expec', 'u', 'y', 'r'||


******************************************************
*                                                    *
* Part 8.  Compute AIC & SIC Minimizing Lag Lengths  *
*                                                    *
******************************************************

*set oshock1 / = oshock(t-1)
*set oshock2 / = oshock(t-2)


source(noecho) c:\myfiles\papers\news\code\aicsic.src


@aicsic(more, nosuppress) sbeg1  send1  4  lagchoice(1,1) lagchoice(1,2)
# pdot expec u y r
*# constant oshock oshock1 oshock2
# constant


*display;display;



                 ***************************************
                 *                                     *
                 *         Begin Bootstrap Routine     *
                 *                                     *
                 ***************************************



*********************************************************************************************
* This program produces the following:
*   IRPOINT_ = NVARS x NSHOCKS RECT[SERIES] of impulse response point estimates,
*             where the response of the ith variable to a shock to the jth equation is
*             IRPOINT(i,j)
*   LBAND_, UBAND_ =
*             NVARS x NSHOCKS RECT[SERIES] of lower and upper points of the confidence
*             interval, each dimensioned in the same way as IRPOINT
* Thus, the impulse response and confidence intervals for y(i) in response to a shock to
* equation j are given by IRPOINT_(i,j) and LBAND_(i,j), UBAND_(i,j).
*
* The (i,j) refer to the order in which the variables are listed in Part 4.a - 4.b (which is
* the same order that you must construct the VAR in Part 6), not the Cholesky ordering defined
* in Part 5.
*
*
*
* Note:  all bootstrap parameters and data types use "_" as the last character of the name
*
*
**********************************************************************************************



**********************************
*                                *
* Bootstap Parameters            *
*                                *
**********************************

compute sbeg1_     = sbeg1+lagchoice(1,LS1); * first period to use to estimate the VAR
                                             * (all lagged values MUST be defined at this date)

compute send1_     = send1;                  * last  period to use to estimate the VAR

compute subtitle_  = %datelabel(sbeg1_)+'-'+%datelabel(send1_)

compute ndraws_    = 1000;   * number of draws

compute nvars_     = nvar+1; * number of variables in the model = VAR eqns + identities
                             * (you must include at least one identity, thus nvars > nshocks)

compute nshocks_   = nvar;      * number of shocks in the VAR (i.e., nvars minus number of identities)
                             * (this is the number of equations in the VAR, not incl. identities)

compute nlags_     = lagchoice(1,LS1);   * number of lags in the VAR

compute nirsteps_  = 25;     * number of impulse response steps to compute

compute width1_     = 0.90;   * enter width for confidence interval 1 (0.95 means 95% intervals)
compute width2_     = 0.80;   * enter width for confidence interval 2 (0.95 means 95% intervals)
compute width3_     = 0.68;   * enter width for confidence interval 3 (0.95 means 95% intervals)


compute [vector] percentiles1_ $
                    = ||(1.0-width1_)/2.0,(1.0+width1_)/2.0||;  * percentiles for boostrapped error bands 1

compute [vector] percentiles2_ $
                    = ||(1.0-width2_)/2.0,(1.0+width2_)/2.0||;  * percentiles for boostrapped error bands 2

compute [vector] percentiles3_ $
                    = ||(1.0-width3_)/2.0,(1.0+width3_)/2.0||;  * percentiles for boostrapped error bands 3


**********************************
*                                *
* Dimension Objects              *
*                                *
**********************************


declare vector[series]        path_(nshocks_) actser_(nvars_)
declare rectangular           pathmatrix_(send1_-sbeg1_+1,nshocks_) irall_(ndraws_,nirsteps_*nvars_) $
                                                                    iroil_(ndraws_,nirsteps_*nvars_)

declare rectangular[series]   ssim_(ndraws_,nvars_) lband1_(nvars_,nshocks_) uband1_(nvars_,nshocks_)  $
                                                    lband2_(nvars_,nshocks_) uband2_(nvars_,nshocks_)  $
                                                    lband3_(nvars_,nshocks_) uband3_(nvars_,nshocks_)  $
                                                    irpointoil_(nvars_,1)                              $
                                                    lbandoil1_(nvars_,1)     ubandoil1_(nvars_,1)      $
                                                    lbandoil2_(nvars_,1)     ubandoil2_(nvars_,1)      $
                                                    lbandoil3_(nvars_,1)     ubandoil3_(nvars_,1)



declare vector[integer]       sbegs_(nvars_)




*****************************************************************
*                                                               *
* Place All Relevant Data Into VEC[SER] = ACTSER_               *
*                                                               *
*   Note 1. The ordering does not affect the Cholesky Decomp    *
*   Note 2. The ordering only affects the position of the VAR   *
*           equations in the model                              *
*   Note 3. You must have nvars entries, VAR equation variables *
*           first, then variables defined by identity           *
*                                                               *
*****************************************************************

* 4.A.  VAR Equation Variables (these are the dependent variables of the VAR)

set actser_(1)    /  = pdot
set actser_(2)    /  = expec
set actser_(3)    /  = u
set actser_(4)    /  = y
set actser_(5)    /  = r


do i = 1, nshocks_
   labels actser_(i)
   # varnames(i)
end do i


* 4.B.  Variables Defined By Identity (these are the dependent variables in the identities)

set spreadr      /  = r(t) - pdot(t)

set actser_(6)    /  = spreadr(t)

* 4.C. Create Vector Of Starting Dates For Each ACTSER_(i)

do i = 1,nvars_
   inquire(series=actser_(i)) sbegs_(i) scrap
end do i


***************************************
*                                     *
* Define Cholesky Ordering            *
*                                     *
***************************************

compute [vector[integer]] corder_ = ||2,3,1,4,5||

compute choleskystring_ = ' / Order = Expec, U, Pdot, y, R'


**********************************************************
*                                                        *
* Define The VAR                                         *
*                                                        *
*    Note 1. List variables in the order given in 4.A    *
*                                                        *
**********************************************************


system(model=VAR1eqns_)
     variables  actser_(1) actser_(2) actser_(3) actser_(4) actser_(5)
     lags 1 to nlags_
     det constant oshock{0 to nlags_} gshock{0 to nlags_}
     *det constant
end(system)



***********************************************
*                                             *
* Define All Identities                       *
*                                             *
*                                             *
* Note: you must hace at least one identity   *
***********************************************



* Interest Rate Spread Implied by the VAR

equation(identity) spreadr_id   actser_(6)
# actser_(5)  actser_(3)
associate          spreadr_id
# 1.0        -1.0



*********************************************************************
*                                                                   *
* Estimate The VAR Model                                            *
*         Capture Residuals                                         *
*         Capture Coefficients                                      *
*         Compute Cholesky Factor                                   *
*         Group Estimated VAR Equations & Identities Into a Model   *
*                                                                   *
*********************************************************************

estimate(residuals=resid_,coeffs=beta_)  sbeg1_  send1_

compute cfactor_ = %psdfactor(%sigma,corder_)


group ID1model_  spreadr_id

compute [model] model1_     = VAR1eqns_ + ID1model_


*******************************************************************************
*                                                                             *
* Compute & Capture Impulse Response Point Estimates, Variance Decomps,       *
* Historical Decompositions & Structural Shocks                               *
*                                                                             *
* Note 1.  The impulse response point estimates are in the NVARS_ x NSHOCKS_  *
*          RECT[SERIES] called IRPOINT_(i,j), i = 1,...,NVARS_ and            *
*          j = 1,...,NSHOCKS_.  The response of variable i to a shock         *
*          to equation j is IRPOINT_(i,j), where the number is that given in  *
*          Part 4A-4B, not by the Cholesky chosen in Part 5                   *
*                                                                             *
*******************************************************************************

display 'Cholesky Factor = ' cfactor_

****************************************
* Compute Structural Errors
*    = S_ERROR, a nshocks_ x 1 VEC[SER]
****************************************

declare vector[integer]  residlist
declare vector[series]   s_error(nshocks_)

enter(varying) residlist
# resid_(1)
do i = 2, nshocks_
   enter(varying) residlist
   # residlist resid_(i)
end do i

make rf_errormat  sbeg1_  send1_
# residlist

compute s_errormat = (rf_errormat)*(tr(inv(cfactor_)))

do i = 1, nshocks_
   set s_error(i) sbeg1_ send1_ = s_errormat(t-(sbeg1_-1),i)
end do i



*******************************************************************
* Compute: Var Decomp, Hist Decomp, Impulse Responses
* Note: HISTORY command will not work on a model with identities
*******************************************************************

errors(model =model1_,decomp=cfactor_)                      *  nirsteps_  *
*history(model=VAR1eqns_,decomp=cfactor_,results=hist_,noadd)  *  (send1_-sbeg1_+1) sbeg1_
impulse(model=model1_,decomp=cfactor_,results=irpoint_)     *  nirsteps_  *

*print(dates) / hist_

******************************************************************************************
*         PATCH:  divide all impulse responses by contemporaneous own response
******************************************************************************************
do j = 1, nshocks_
   compute rescale = irpoint_(j,j)(1)
   do i = 1, nvars_
      set irpoint_(i,j) 1 nirsteps_ = irpoint_(i,j)/rescale
   end do i
end do j
******************************************************************************************



*******************************************************************************
*                                                                             *
* Compute & Capture Impulse Response Point Estimates to Oil Dummy Shock       *
*                                                                             *
* Note:  we compute this as the difference between 2 forecasts, one of which  *
*        has a (temporary) unit increase in Hamilton's oil dummy measure in   *
*        the first period of the forecast.                                    *
*                                                                             *
*******************************************************************************

****************************
* Compute Baseline Forecast
****************************

forecast(model=model1_,results = noshockfor_) * nirsteps_  sbeg1_

***********************************************
* Compute Forecast With Higher Dummy Variable
***********************************************

set oshock sbeg1_ sbeg1_ = oshock + 1.0

forecast(model=model1_,results = shockfor_)   * nirsteps_  sbeg1_

set oshock sbeg1_ sbeg1_ = oshockreplace


*************************************
* Compute Implied Impulse Responses (as the difference between the forecasts)
*************************************

do i = 1, nvars_

   set irpointoil_(i,1) / = shockfor_(i) - noshockfor_(i)

end do i

clear shockfor_  noshockfor_



*******************************************************************************
*                                                                             *
* Compute Simulated VAR Data                                                  *
*                                                                             *
* Note 1. The simulated values are in the NDRAWS x NVARS RECT[SERIES]         *
*         called SSIM_(draw,i),where the row (draw) indicates the replication *
*         number and the column (i) indicates the variable simulated          *
*                                                                             *
*******************************************************************************

* Conduct Simulation

smpl sbeg1_ send1_
do draw_ = 1, ndraws_

   boot entries_ / sbeg1_ send1_

   do i = 1,nshocks_
      set path_(i) / = resid_(i)(entries_(t))
   end do i

   do i = 1, nshocks_
      do j = 1, send1_-sbeg1_+1
         compute pathmatrix_(j,i) = path_(i)(sbeg1_-1+j)
      end do j
   end do i

   forecast(model=model1_,matrix=pathmatrix_, results=results_)

   do i = 1,nvars_
      set ssim_(draw_,i) / = results_(i)
   end do i

end do draw_


smpl
* Reset SMPL Above This Line


*********************************************************************************************
*                                                                                           *
* Estimate VAR  on Simulated Data                                                           *
*          Compute Impulse Responses                                                        *
*          Calculate CIs                                                                    *
*                                                                                           *
* Note 1.  The lower and upper confidence intervals are given in the NVARS x NSHOCKS        *
*          RECT[SERIES] called LBAND_(var,shock) and UBAND_(var,shock), where the row (var) *
*          indicates the response variable and the column (shock) indicates the equation    *
*          that is shocked. Note that the variable number and equation-shock number are     *
*          given in Part 4.A-4.B, not the Cholesky ordering in Part 5.                      *
*                                                                                           *
*********************************************************************************************

* Splice Initial Conditions (given by history) On Top Of Simulated Values

do draw_ = 1, ndraws_
   do var_ = 1, nvars_

      set ssim_(draw_,var_) sbegs_(var_) sbeg1_-1  = actser_(var_)

   end do var_
end do draw_


* Impulse Responses

do shock_ = 1,nshocks_

do draw_=1,ndraws_

   * Estimation

   system(model=VAR2eqns_)
     variables ssim_(draw_,1) ssim_(draw_,2) ssim_(draw_,3) ssim_(draw_,4) ssim_(draw_,5)
     lags 1 to nlags_
     det constant oshock{0 to nlags_} gshock{0 to nlags_}
     *det constant
   end(system)

   estimate(noprint,noftests,coeffs=simbeta_) sbeg1_ send1_



   * Compute Coefficient Averages (over draws)[Not Used In This Program]

   if draw_.eq.1
     {
       compute musimbeta_ = (1.0/ndraws_)*simbeta_
     }
   else
     {
       compute musimbeta_ = musimbeta_ + (1.0/ndraws_)*simbeta_
     }

   * Define Identities

   * Interest Rate Spread
   equation(identity) spreadr_id   ssim_(draw_,6)
   # ssim_(draw_,5)  ssim_(draw_,3)
   associate          spreadr_id
   # 1.0           -1.0


   group ID2model_  spreadr_id

   compute [model] model2_     = VAR2eqns_ + ID2model_


  * Compute Impulse Responses

   * Note 1: the equation shocked by the loop index, called "shock", is given by the order of
   *         equations in the VAR, not the order given by the cholesky factorization
   * Note 2: IR is an Nvars x 1 RECTANG[SER], with the ordering given by the ordering of the
   *         equations in the VAR, not the order of the cholesky factorization, and the
   *         ordering of the identities


   compute cfactor_ = %psdfactor(%sigma,corder_)
   impulse(model=model2_,noprint,decomp=cfactor_,results=ir_) * nirsteps_ shock_




   *****************************************************************************************
   *        PATCH:  divide all impulse responses by contemporaneous own response
   *****************************************************************************************
   compute rescale = ir_(shock_,1)(1)
   do i = 1, nvars_
     set ir_(i,1) 1 nirsteps_ = ir_(i,1)/rescale
   end do i
   *****************************************************************************************




   * Capture Impulse Responses For This Draw


   do var_ = 1,nvars_
      do step_ = 1,nirsteps_
          compute irall_(draw_,step_+(var_-1)*nirsteps_) = ir_(var_,1)(step_)
      end do step_
   end do var_

  * Compute & Capture Impulse Responses For Oil Dummy

   if shock_.eq.1
   {

      * Compute Baseline Forecast

      forecast(model=model2_,results = noshockfor_) * nirsteps_  sbeg1_

      * Compute Forecast With Higher Dummy Variable

      set oshock sbeg1_ sbeg1_ = oshock + 1.0

      forecast(model=model2_,results = shockfor_)   * nirsteps_  sbeg1_

      set oshock sbeg1_ sbeg1_ = oshockreplace

      * Compute Impulse Response As Difference Between The Forecasts & Capture

      do var_ = 1,nvars_
         do step_ = 1,nirsteps_
            compute iroil_(draw_,step_+(var_-1)*nirsteps_) = shockfor_(var_)(sbeg1_-1+step_)   $
                                                           - noshockfor_(var_)(sbeg1_-1+step_)
         end do step_
      end do var_

      clear noshockfor_ shockfor_

   }

end do draw_


   * Compute Confidence Intervals For This VAR Shock & The Oil Dummy Shock (by variable and step)

   do var_ = 1,nvars_
      do step_ = 1,nirsteps_

         * VAR Shock Responses
         compute col_ = %xcol(irall_,step_+(var_-1)*nirsteps_)
         compute lowerupper1_ = %fractiles(col_,percentiles1_)
         compute lowerupper2_ = %fractiles(col_,percentiles2_)
         compute lowerupper3_ = %fractiles(col_,percentiles3_)


         set lband1_(var_,shock_) step_ step_ = lowerupper1_(1)
         set lband2_(var_,shock_) step_ step_ = lowerupper2_(1)
         set lband3_(var_,shock_) step_ step_ = lowerupper3_(1)


         set uband1_(var_,shock_) step_ step_ = lowerupper1_(2)
         set uband2_(var_,shock_) step_ step_ = lowerupper2_(2)
         set uband3_(var_,shock_) step_ step_ = lowerupper3_(2)

      end do step_
   end do var_

end do shock_



*********************************************
*                                           *
* Graphs of Impulse Responses:  VAR Shocks  *
*                                           *
*********************************************


*********************
* Expectations Shock
*********************

spgraph(vfields=4,hfields=2,header='Expectations Shock')

   graph(header='Inflation Response',nodates,number=0) 7
   # lband1_(1,2)     1 nirsteps_   1
   # irpoint_(1,2)    1 nirsteps_   4
   # uband1_(1,2)     1 nirsteps_   1
   # lband2_(1,2)     1 nirsteps_   3
   # uband2_(1,2)     1 nirsteps_   3
   # lband3_(1,2)     1 nirsteps_   6
   # uband3_(1,2)     1 nirsteps_   6

   graph(header='Expectations Response',nodates,number=0) 7
   # lband1_(2,2)     1 nirsteps_   1
   # irpoint_(2,2)    1 nirsteps_   4
   # uband1_(2,2)     1 nirsteps_   1
   # lband2_(2,2)     1 nirsteps_   3
   # uband2_(2,2)     1 nirsteps_   3
   # lband3_(2,2)     1 nirsteps_   6
   # uband3_(2,2)     1 nirsteps_   6

   graph(header='Unemployment Response',nodates,number=0) 7
   # lband1_(3,2)     1 nirsteps_   1
   # irpoint_(3,2)    1 nirsteps_   4
   # uband1_(3,2)     1 nirsteps_   1
   # lband2_(3,2)     1 nirsteps_   3
   # uband2_(3,2)     1 nirsteps_   3
   # lband3_(3,2)     1 nirsteps_   6
   # uband3_(3,2)     1 nirsteps_   6

   graph(header='Output Response',nodates,number=0) 7
   # lband1_(4,2)     1 nirsteps_   1
   # irpoint_(4,2)    1 nirsteps_   4
   # uband1_(4,2)     1 nirsteps_   1
   # lband2_(4,2)     1 nirsteps_   3
   # uband2_(4,2)     1 nirsteps_   3
   # lband3_(4,2)     1 nirsteps_   6
   # uband3_(4,2)     1 nirsteps_   6

   graph(header='Interest Rate Response',nodates,number=0) 7
   # lband1_(5,2)     1 nirsteps_   1
   # irpoint_(5,2)    1 nirsteps_   4
   # uband1_(5,2)     1 nirsteps_   1
   # lband2_(5,2)     1 nirsteps_   3
   # uband2_(5,2)     1 nirsteps_   3
   # lband3_(5,2)     1 nirsteps_   6
   # uband3_(5,2)     1 nirsteps_   6

spgraph(done)



***********************************
*                                 *
* Graphs of Structural Shocks     *
*                                 *
***********************************


*spgraph(vfields=3,hfields=2,header='Estimates of Structural Shocks', $
*        subheader=subtitle_+choleskystring_)

*    graph(header='Inflation Shocks')       1
*    # s_error(1)   sbeg1_  send1_ 1

*    graph(header='Expectations Shocks')    1
*    # s_error(2)   sbeg1_  send1_ 1

*    graph(header='Unemployment Shocks')    1
*    # s_error(3)   sbeg1_  send1_ 1

*    graph(header='Monetary Policy Shocks') 1
*    # s_error(4)   sbeg1_  send1_ 1

*    graph(header='Commodity Price Shocks') 1
*    # s_error(5)   sbeg1_  send1_ 1

*spgraph(done)


**************************************************************
*                                                            *
* Output Impulse Responses & Hist Decomps to a Worksheet     *
*                                                            *
**************************************************************

*******************************
* VAR Shock Impulse Responses
*******************************

set steps 1 nirsteps_ = t-1

* Confidence Intervals 1

open copy c:\myfiles\papers\news\code\output\March08\ordered_first\new_timing\6months_ahead\dec09\revisions\VarIR_LS5_BaseWDUMMY_90_gdprt_12mo.xls
copy(nodates,format=xls,org=obs) 1  nirsteps_  $
                                 steps  irpoint_  lband1_   uband1_

* Confidence Intevals 2

*open copy d:\myfiles\papers\sunspots\Tom's code\output\VarIR_PreV_Prg7B1_BaseO_80.xls
*copy(nodates,format=xls,org=obs) 1  nirsteps_  $
*                                 steps  irpoint_  lband2_   uband2_

* Confidence Intervals 3

open copy c:\myfiles\papers\news\code\output\March08\ordered_first\new_timing\6months_ahead\dec09\revisions\VarIR_LS5_BaseWDUMMY_68_gdprt_12mo.xls
copy(nodates,format=xls,org=obs) 1  nirsteps_  $
                                 steps  irpoint_  lband3_   uband3_


***********************
* Historical Decomps
***********************

*open copy c:\myfiles\papers\sunspots\Tom's code\output\Hist_PreV_Prg7B1_BaseO.xls
*copy(dates,format=xls,org=obs) /                           $
*                                  pdot  epdot  u  r  pcom  $
*                                  hist_





































































































































































































































































































































































































